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Abstract. 

In order to elucidate the relationship between rate-independent hysteresis and 
metastability in disordered systems driven by an external field, we study the Gaussian 
RFIM at T = on regular random graphs (Bethe lattice) of finite connectivity z and 
compute to 0(1/ z) (i.e. beyond mean-field) the quenched complexity associated with 
the one-spin-flip stable states with magnetization m as a function of the magnetic field 

H. When the saturation hysteresis loop is smooth in the thermodynamic limit, we 
find that it coincides with the envelope of the typical metastable states (the quenched 
complexity vanishes exactly along the loop and is positive everywhere inside). On 
the other hand, the occurence of a jump discontinuity in the loop (associated with an 
infinite avalanche) can be traced back to the existence of a gap in the magnetization of 
the metastable states for a range of applied field, and the envelope of the typical 
metastable states is then reentrant. These findings confirm and complete earlier 
analytical and numerical studies. 
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1. Introduction 

At low temperature, disordered systems exhibit a large number of metastable states 
which are responsible for their slow relaxation dynamics and their irreversible, hysteretic 
response to an external field. In many cases this response is very well reproducible 
and independent of the rate at which the field is changed, provided it is low enough. 
This shows that these systems always explore the same set of metastable states when 
driven by the external field and that thermally activated processes are irrelevant on 
experimental time scales. Classical examples of such behavior are magnetic materials. 
When the external magnetic field is slowly cycled from large negative to large positive 
values and back, the induced magnetization displays a (saturation) hysteresis loop with 
a characteristic shape that depends on the mechanisms for the formation and evolution 
of the magnetic domains[lJ. Other field histories give rise to a complicated subloop 
structure that is also well reproducible. 

By definition, the future evolution of a hysteretic system depends on its past 
history, and the theory of rate-independent hysteresis generally focuses on a dynamical 
description in order to describe history-dependent effects such as return-point memory [TJ 
El [3]. However, the hysteresis loop may also be viewed as a special path among the 
metastable states and one may wonder whether it is possible to characterize this path 
without having to follow the dynamical evolution. Surprisingly, a general answer to 
this question is not yet available. In the case of attractive (ferromagnetic) interactions, 
however, an important fact is known: if the dynamics obeys the so-called no-passing 
rule [2] (i.e. preserves some natural partial ordering of the microscopic configurations), 
there are no metastable states outside the saturation hysteresis loop[l]. Since, as a 
general rule, the number of metastable states scales exponentially with the system size, 
this result implies that the corresponding entropy density (or complexity) is not positive 
outside the loop. What still remains to be understood is the situation inside the loop. 

A convenient framework for studying this problem is the zero-temperature Random 
Field Ising model (RFIM) driven by an external field. This is a model for disordered 
ferromagnets which is relevant to a large class of materials where disorder induces 
random fields [5J. With the standard Glauber dynamics (which does satisfy the no- 
passing rule at T = 0[2]), the metastable states are the so-called single-spin-flip stable 
states, each state being characterized by its energy and magnetization. The quantity 
to be studied in relation to hysteresis is Af(m,H), the number of metastable states at 
the field H with given magnetization per spin m. One can then reformulate the above 
question as follows: in which region of the field-magnetization plane is the logarithm 
of J\f(m, H) an extensive quantity ? In particular, does the saturation hysteresis loop 
identifies with the boundary of this region ? To answer these questions, one must study 
the properties of a typical system and compute the quenched complexity £<g(m, H) of 
the metastable states[l]. This is a nontrivial analytical or numerical task in general, 
and there exists very few results on the typical number of metastable states in finite- 
connectivity disordered spin systems. For the Gaussian RFIM, the only analytical 
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result so far has been obtained in one dimension [1] (or, equivalently, on a Bethe lattice 
of connectivity z — 2), showing that the curve Sg(m, H) = coincides with the 
hysteresis loop. Numerical results on the Bethe lattice of coordination z = 4 and on 
the cubic lattice [6] strongly suggest that this property always holds when the hysteresis 
loop is smooth in the thermodynamic limit (which occurs when the disorder is large 
enough). This means that the typical number of metastable states is exponentially 
large everywhere inside the loop. The situation is more complicated when the loop has 
a jump at some critical value of the field, which corresponds to an out-of-equilibrium 
phase transition associated with the appearance of a macroscopic "avalanche" [2] . This 
happens at low disorder in three and higher dimensions^ E] or on a Bethe lattice for 
z > 4[8|. It has been suggested[l] that this jump reveals the existence of a gap in the 
magnetization of the typical metastable states in a finite range of H, gap which may also 
explain the reentrant hysteresis loops observed when changing the driving mechanism [9 J. 
This assumption is also supported by numerical computations [6], but again there is no 
analytical proof. The aim of the present work is to give further evidence to this scenario 
(both above and below the critical disorder) by computing analytically the quenched 
complexity of the Gaussian RFIM on a Bethe lattice to order 1/z, and comparing with 
the exact expression of the hysteresis loop [8] at the same order. This will also help us 
to clarify the approach to the mean-field fully-connected limit |10j. 

The paper is organized as follows. In section II, we define the model and present 
the general formalism for calculating the quenched complexity. Section III is devoted to 
the mean- field limit which is recovered when z — > oo. The calculation of the complexity 
at the order 1/z is presented in section IV and the numerical results are discussed in 
section V. We finish with a brief discussion and some perspectives in section VI. More 
details on the analytical calculations are given in the appendices. 

2. General formalism 

We start with the RFIM Hamiltonian 



describing a collection of N Ising spins (sj = ±1) placed on the vertices of a graph 
of connectivity z and submitted to an external uniform field H. A ferromagnetic 
interaction J > couples the spin i to its z neighbours j and {hi} is a collection 
of random fields drawn identically and independently from a Gaussian probability 
distribution V(h) with zero mean and standard deviation A. In the zero-temperature 
system with metastable dynamics[2], the field H is adiabatically varied (for instance 
from — oo to +oo and back) and each spin is aligned with its local effective field 




(1) 



<i,j> 




(2) 



so that 




(3) 
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In the thermodynamic limit, the resulting hysteresis loop m(H) (where m = (1/iV) £\ Si 
is the magnetization per spin) is smooth for A > A c (z) and discontinuous for A < A c (z), 
where A c (z) is a critical value of the disorder at which avalanches of all sizes occur 
(an avalanche corresponds to a jump in the magnetization curve m(H)). A complete 
analytical description of this behavior is available in the mean-field model[2], [7] which 
can be obtained for instance by placing the spins on a fully-connected lattice. In this 
limit, there is no hysteresis above A° = lim^oo A c (z) = Jy^2/n and the number of 
metastable states is not exponentially large [TU] (but it may be larger than 1, even along 
the so-called "unstable" branch for A < A°). 

Our goal in this paper is to compute the first non-zero term in the 1/z expansion 
of the complexity Sg(m, H). Although the corrections to the mean-field result are 
certainly different on the Bethe lattice and on the hypercubic lattice, we only consider 
the former case since an analytical expression of the hysteresis curve is available for 
any z[5]. This also allows us to use the machinery developped in Ref. [J] to compute 
Sg(m, H). The possible influence of geometry will be briefly discussed in the conclusion. 

Sg(m, H) is defined as 

1 , 1 



E Q (m,H) = lim — lnJVYm, H) = lim - lim — M(m,H) n - 1 (4) 

where, as usual, the order of the limits iV — > oo and the number of replicas n — > 
has been inverted. The starting point of our calculation is the following expression of 
Af(m, H) n in the large- iV limit [I]: 

/n 
\\ dg a H dc{(T, T )e N ^ {chg) (5) 
1 /T T 



a=l <T,T 



where 



and 



it 



Hie}, 9) = A({c}, 0) + §-f£ <a, r)c(r, cr) - m£ g a (6) 



2 2 

err a=i 



A({c},g) = \nJ2 [ V{h)dh j d*dy e iy ^- n - h) [ ^ c(r, cr) e - ay - T ] * f[ e sV °e(iV a ) 

CT T a=l 

(7) 

where Q(x) is the Heaviside function. In these equations, the bold symbols denote 
n-component vectors in replica space, e.g. x = {x a ; a = l,...n}, (with h = hi and 
H = HI); c(cr, r) is an order parameter, function of the two binary vectors a and 
r [o a = ±1, r a = ±1), which allows the decoupling of the sums over the sites i. It 
generalizes the single-argument order-parameter function c(cr) introduced in Ref. [TTJ 
to study finite-connectivity spin glasses (see also Ref.[T2]). g = {g a } is a set of Lagrange 
multipliers associated with the constraint ^\ s" = Nm in each replica. (The interested 
reader is referred to Ref. [1] for the derivation of these equations [13J. Apart from these 
three equations, the present paper is self-contained.) 
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In the large- iV limit, the integral in Eq. (j^D can be evaluated by the method of 
steepest descent. The 2 2n order parameters c(cr, r) and the n Lagrange mutipliers g a 
are determined through the saddle-point equations that correspond to the extremization 
of the functional F({c},g), 

c(<r , T) .I^£M (8) 

z oc(t, cr) 

_^|M, (9) 

where a is any of the replicas 1, ...,n. One can easily see that the solution of Eq. ([S]) 
(that will be denoted as c*(cr,r), dropping the dependence on g and H) satisfies the 
normalization condition 

^c*(<x,t)c*(t,<t) = 1 (10) 

(T.T 

so that Eqs. ([9]) can be rewritten as[l] 

m = ^a a c>,T) C *(T,<r) . (11) 

(T.T 

Assuming that the solution of Eq. (TlTI) in the limit n — > is the same for all 
replicas, namely g* a = g*(m, H), we expand A*(m, H) = A(c*(cr, t), p*) in powers of n, 

A*{m,H) = A<°\m,H) + nA< 1 \m,H) + 0(n 2 ) , (12) 

which, after insertion in Eq.(jSJ), yields 

E Q (m, H) = A* (1) (m, H) - mg*(m, H) , (13) 

where we have used that A*^(m, H) is equal to zero, which is necessary to obtain a 
well-defined n — >• limit. 

As discussed in Refs. jU [6] , it is in fact convenient to treat the Lagrange multiplier 
g as a free parameter and to consider that the magnetization is a function of g given 
by the solution of Eq. ffTTT) in the limit n — » 0. Defining A*(g,H) = A(c* (cr , t) , g) , we 
then write 

£ Q (<7, if) = £ Q (m*(<7, H),H) = A«%, H) - m*(g, H)g (14) 

where m*(g, H) = lim n ^ J2cr r cr<lc *( cr ? t )c*(t, cr) = <9A* (1 )(g, H)/dg. In other words, 
Eg(m, .ff) and A*^\g,H) are mutually connected by a Legendre transform with 
—g*(m,H) being the slope of the complexity curve, i.e. g*(m,H) = —dT,Q(m,H)/dm. 
It is crucial that the mapping {m*(g, H),T,q(g, H)} \— > Sg(m, H) is unambiguously 
defined, as will be discussed in section V. Note also that the total complexity at 
the field H, irrespective of the magnetization, is obtained by setting g = in 
Eq. ffMl) . This corresponds to the maximum of Eg(m, H) that occurs at a magnetization 
m max (H) = m*(g = 0,H). 



The T = random-field Ising model on a Bethe lattice with large coordination number: hysteresis and metasi 

We now turn to the 1/z expansion of Eq(?7j, H). Using the rescaling Jz — > J, we 
rewrite Eq. (jTJ) as 



A({c},g) =lnJ2 e9S J V{h)dh J dxdy e^ x ' y) Y[e(x a o- a ) 



(15) 



where s = £ a a a and 

0(x, y) = iy.(x - H - h) + z In ^ c(r, <r) e -^ y - T . (16) 

T 

We then formally expand Eqs. (jSJ) and fflOl) in powers of 1/z, assuming that 

A*(<?, #) = zA^g, H) + A*G?, H) + -Afo, if) + ... (17) 

z 

m*(g, H) = m* (g, H) + -ml{g, H) + ... (18) 



and 



1 1 

c*(tr, t) = c*(a, t) + -c*(<r, r) + -^(<r, r) + ... (19) 



Z ' z 2 



This yields a set of coupled equations to be solved at each order in 1/z. The main 
difficulty of the calculation is that successive orders are not independent. Indeed, one 
has 

0(x, y) = z In ( c* {t, trj) + o (x, y) + -^(x, y) + -^ 2 (x, y) + ...(20) 

T 

with 

0o(x, y) = iy.(x — H — h) + lK ^ J y oK 1 n (21a) 

s = E T ic*2(T,<T) -iJy.TcKr,*) - f (y .t) 2 c*(t , a)] 1 r ^tcKT.^-ay.TC^T,^] ^ 
<Pl[ * ,y) Er^a) 2l £ T c8(T,<r) J 

(216) 

etc.... Inserting Eq. ffTBI) in Eq. (jSJ) and using these expansions, we obtain 

A*_ 1 + ln^c*(r, CT ) = (22) 

T 

at leading order, and 

c*(<r, r) = e ~ A ° +9S J V{h)dh J dxdy e* o(x ' y) JJ 6(xV a ) , (23) 

a 

c *( ct ,t) = e - A *> +9S / V(h)dh J dxdy e* o(x ' y) JJ 6(a;V a ) [ - A* + 0i(x,y) - iJy.r 

J J a 

-J2[ct(T\a)-Uy.T> C *(T>,*)] (24) 

T' 
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at the two next orders. Similarly, the expansion of the normalization equation, Eq. ([TUj) . 
yields 

J>*(<t,t)c*(t,<t) = 1, (25) 

(T.T 

J2<%(<7,t)c* 1 (t,<t) = 0, (26) 

(T T 

etc... (Note that the dependence on g and/or H is not always indicated in order to 
simplify the notations. We shall do the same in the following, except when needed 
explicitly.) 



3. Mean-field limit (z — > oo) 

We first consider the limit z — > oo and check that the results of the mean-field 
model[2j [7J [10] are recovered. The calculation, which in this setting is not trivial, 
proceeds as follows. Inserting Eq. (21a) into Eq. (1231) and performing the integrations 
over y and x, we obtain 



nl „. , , r -K+9s+T, T icl(T',(T)/j: Tl c*(T',cr) 



Cn C, T 



dhV(h)Y[e(x* a a a ) (27) 



with 



Since the r.h.s. of Eq. (1271) does not depend on r, Cg is a function of its first argument 
only that will be denoted as Cq(ct). As a result, the quantity x* a does not depend on cr. 
It also does not depend on a if one assumes that all replicas are equivalent. Therefore, 
because of the integration over h, the only functions that differ from zero are Cq(1) and 
Cq(— 1) (i.e. all a a are equal to 1 or —1). From Eq. (1231) . these two functions satisfy 

cg(l) + cg(-l) = 1 (29) 

(the other solution of Eq. (125]) . Cq(1) + Cq(— 1) = —1, is not physically relevant). Eq. 
( f29l) is actually a prerequisite for having a finite limit when z — >• oo, since it implies 
from Eq. (122!) that A* j = in the expansion (TTTj) . Defining Xi(er) = ^ r c^(r, cr) and 
Ac* = cj(l) - cg(-l), we then rewrite Eq. ([271) and (E2} as 

i[l ± Ac*] = e -AS±s»H-Xi(±l) J dhV(h)Q(±x*) (30) 



and 



where 



^^toJWHWHe^-VjWKM-S) (3!) 
x * = H + h + JAc* . (32) 
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This is not yet a closed set of equations for computing Ac*, and Aq because the two 
quantities Xi(±l) depend on c*(<r, r). Therefore, the equations at the next order in 
1/z, Eqs. ( 1241 and ( 1261) . must be also considered. The latter can be written as 

(1 + Ac*)Xx(l) + (1 - Ac5)Xi(-l) = (33) 

so that we only need to extract from Eq. ( [241 the quantity AXi = X\(l) — X\{— 1). A 
quick look at this equation shows that the solution takes the form 

cft<r, r) = F(a) + e ~ A o + v s+x ^I((T, r) (34) 

where 



) (35) 



y 



I((T,r) = -iJ j V{h)dh I dxdy e iy ' M,) y.r[]e(xV a 

J J a 

and F(cr) groups all other terms. The simple result follows: 
AXi = ^e- A *' +£ ' s+Xl ( cr )[/(cr, 1) - I((T, -1)] 

= -2ijJ2 e ~ A * 0+9S+Xl(CT) I P(h)dh j rfxJj0(xV a ) f dy e iy - (x - x * } 

<7 J J a J a 

(36) 

Using —i^2 a y a = 9e 4y ^ x ~ x *- ) /dx* = <9e ty '( x ~ x ** ) /dh and integrating by parts over h, we 
find that all <r a 's must be equal in the sum over cr, which yields 

AX 1 = 2JV*e- A ° [ e 9 n+x ^ - e -s«+^i(-l)] (37) 

where V* is a shorthand for V(-H - JAc* ). Eqs. (J30D, ([31]), ([33]), and (1371) 
now form a closed system and Aq, Ac*,, and Xi(±l) can be calculated in the limit 
n —y by expanding in powers of n (i.e. assuming that A* = A W + nA^> + 
Ac* = Ac* {0) +nAc* {1] + ..., and X 1 (±l)=x{ 0) (±l)+nXj 1) (±l)...). At the lowest order 
in n, this gives 

A<; (0) = (38) 
X; o) (±l)=0 (39) 
Ac {0) = 2p{Acf ] ) - 1 (40) 

where p(m) = f™ H _ Jm V(h)dh = (1/2)[1 + erf([# + Jmj/Av 7 ^)]. Note that, as stated 
before, A W = is a prerequisite for a proper limit n — - > 0. In addition, Eq. lull readily 
gives 

m (H) = lim *°<Sfo t) c o(t, <t) = Ac* (0) , (41) 



(XT 



so that m*> satisfies the self-consistent mean-field equation [21 171 [TO] 



m, 



\(H) = exf([H + JmS(^)]/A^) (42) 



and does not depend of g. 

At the order n, the solution of Eqs. (13^ . (131]) . (1331. and ([37]) yields A* Q {1 \g,H) = 
gm^H) (so that m*, = dA^ /dg, as it must be), and the corresponding complexity 
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TjQfi(g,H) = AQ^\g,H) — m*Q(H)g is thus identically zero. In the H — m plane, this 
means that the quenched complexity is zero when m = ru^H), i.e. along the mean- 
field loop, and is not defined otherwise. On the other hand, a trivial calculation shows 
that the annealed complexityjl] XU(m, H) = liniAr^oo In Af(m, H) is also zero when 
m = itlq{H) and strictly negative otherwise. We stress that a zero complexity does 
not imply that there is a unique stable state in the thermodynamic limit. The actual 
number of metastable states along the curve m^H) can be larger than 1, even along 
the so-called "unstable" branch for A < A°, as shown in Ref.|10j. 

The expressions of Ac {1) and x[ 1] (±l), which are needed in the subsequent 
calculations, are given at the end of Appendix A. 



4. l/z correction to the quenched complexity 

We now compute the l/z correction to the complexity which, from Eq. (fT4l) . is given 
by 

X Q ,i(g,H) = Af ) (g,H)-ml(g,H)g . (43) 

To compute A*^(g, H) (and then m\{g,H) by derivation), we have to fully solve Eq. 
( |24l) so to obtain the explicit expressions of the functions J(cr, r) and F(cr) defined by 
Eq. (1341 . After straightforward but lengthy algebraic manipulations (which are detailed 
in Appendix A), we find that A* is given by 

At = -~[Xi(l) 2 co(l) + X 1 (-l) 2 c (-l)\ - + Xi(-l)] - -[^(1) + ^(-l)] 1 

-^JV*[Xx(l) - X 1 (-l)][e- A S +9ri+Xl(1) - e -AS-fl*+*i(-l)] 

_1 j2p*2^-A ( *+ 9 n+Xi(l) + e -A*-gn+Xi(-l)j2 

-JV*Ac* {[l + A 1 (l)]e~ A " +9 ' 1+Xl(1) - [1 + A 1 (-l)]e- A S- 9 " +Xl( - 1) } 

+— (1 - Acjf )P'*[e~ A( * +9 ™ +Xl( } - e - A o-9 n+Xl (- 1 )] + - y^[2J~P*e~ A Q +gs+Xl(cr) ] 2 



(44) 



and that X\{tr) satisfies the equation 

Xi(o-) - Xi(l) = 2JP*e- A S[e 9S+Xl(cr) - e 9n+Xl{1) ] . (45) 

This equation plays a central role in our study and will be discussed in detail below. 
Let us first define the quantity T(g,H) = Y J(T [2JV*e-^ +9S+x ^ (T) } 2 = E<r[ x i(<x) - 
Xi(l) + 2JV*e- A ° +9n+Xl ( 1) ] 2 and expand both sides of Eq. flHD in powers of n (with 
T(g, H) = T<®(g, H) + nT^(g, H) + ...). Using Eqs. (40), (41), and flA~T5l) . this yields 

Af\g 1 H) = -2J 2 Vf+ 1 -T^(g) (46) 

and 

A;«( 9 , H) = 9 (1 - 0^(1 - 4 T -^) + ir«( ff ) (47) 
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where P * is a shorthand for V(H + Jm*(H)) (i.e. V* = P * + 0(n)), and V'(h) = 
dV(h)/dh (for brevity, the dependence of T(g, H) on the field H is dropped). From this 
we obtain 

ml(g,H) = dAf\g,H)/dg 

-n m* 2 ) J2V ° n 1 1 I 1 9T(1)( ^ f4tt 

- (1 " m ° ) T^2Ji^ (1 " 4 T^2J^ ) + 2^^" (48) 

and 

^ Q A9,H) = \[T^\g)-g^l]. (49 ) 

(It can be checked through lengthy calculations that the same expression of m\(g : H) is 
recovered by expanding directly Eq. (TTTT) ) to 0(1/ z).) 

The only remaining task is to solve Eq. (H5|) for X±(a) and to calculate T^(g) and 
T^(g). This equation can be conveniently rewritten as 

W(a)e w{(T) = -A(g)e- A{9)+9( - s ~ n) (50) 

where W(o) = -X x (a) + X x (l) - and = 2JV*e- k ^ +gn+Xl{ - 1 \ Because of 
replica symmetry, W(tr) only depends on s = ^ a 0" a and is a function of the single 
variable z(s,g) = -A(g)e~ A ^ +9{ - s ~ n \ Eq. ([50]) becomes 

W(z)e w{z) = z , (51) 

which tells us that W(z) is the so-called "Lambert function" |T4] (remarkably, this 
function already appears in the calculation of the average number of metastable states 
along the mean- field curve rriQ(H)[^). The whole problem thus amounts to compute 

T(g) = J2w 2 (z(s,g)) (52) 

as n — > 0. 

Let us first consider the case g = which, as noted ealier, corresponds to the 
maximum of the complexity and thus yields the total number of metastable states at 
the field H. According to Eq. (|48|) . in order to obtain the corresponding magnetization 
m max (H) = m(g = 0,H), one also needs to compute lim^o dT^(g)/dg. We then 
expand both sides of Eq. floTj) in powers of g, using the formula for the n-th derivative 
of the Lambert function |14j. This allows us to perform the sum over cr in each term of 
the expansion (computing Ylcr s > Scr s2 > etc.). The result is then expanded in powers 
of n, which gives 

T (0 \g) = t\ (53) 



and 



m , x 2l 2t 2 2t 2 t 2 (2 + t) 2 t 2 (2t 2 + 9t- 

T^(g) = t 2 \n2 + L_ '-—g - —g + f ) ■ J g 2 + 1 g/1 , — 7 

(54) 



' 4 

1 + t* l + t* ' (l + t) 3 " ' 6(1 + 1) 7 
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where A(g) = A® + nA^g + ... and t = W{-A^e- A{0) ). Since A^ (0) = x[°\l) = 
(Eqs. ( 1381) and (39)), one has 

A (o) = 2 jp* ; (55) 
and, using Eqs. (IA. 15[) . 

9 TV* T 2 T>'* 
A® = (1 - ml) '° + 2(1 - mf) , 7 - . (56) 

A question now needs to be addressed: which branch of W(,z) should be chosen? 
Recall that the Lambert function has two real branches, Wq{z) and W-i(z), with a 
branch point at z = — l/e[14j. The principal branch Wq(z) takes on values between —1 
to +oo for z > — 1/e (with Wo(0) = 0) and W-i(z) takes on values between — oo and 
— 1 for — 1/e < z < 0. Accordingly one has 

W (-xe~ x ) = -x if x < 1 

W Q {-xe- x ) = -x' if x > 1 (57) 
where x' is the solution (smaller than 1) of the equation x'e~ x ' = xe~ x . Similarly, 
W_i(-xe~ x ) = -x if x > 1 

VVLiC-xe -1 ) = -x' ifx<l (58) 

where x' > 1. 

The choice of the branch is determined by the fact that A* 1 - ) must vanish at all 
orders in 1/z in order to get a proper limit when n — > 0. Eq. (H6l) then implies that 

2^%) = 4J 2 P* 2 = A^ 2 (59) 

whence t = W(— A^e^ Am ) = — A^ in Eq. ( 1531) . Therefore, according to Eqs. ( 1571) and 
(|58|) . one has to choose the principal branch W / o(- 2; ) for < 1 and the branch W-i{z) 
for A^ > 1 (note that Wo is the only branch that comes into play in the calculations 
of Ref . [TO] ) . The series (|54p then becomes 



TW(g) = A^ In 2 + 2A^A^g - -^—g 



2A (o)2 A^ 2 (2-A^) 



AW* (1-A(°)) 3 



2 



^4(o)2 (2A (o)2 _ Q A (Q) _ g) ^ ^ 



6(1 -A(°))< 
so that, from Eqs. (SSJ) and (1491 . 

4 72-p*2 /2-p*' 

m™ ax (H) = -m* + (1 - m* 2 ) , (61) 

and 

T%?{H) = E Qjl (g = 0,H) = 2J 2 P* 2 ln2 . (62) 

This result will be commented on in the next section. More generally, the problem 
is that the series ( l60i) is strongly divergent when g is large. We therefore need to find 
an integral representation that allows the sum to make sense for an arbitrary value of 
g. With this in mind, we write z(s,g) = — A(g)e~ A ^ [1 + (e 9 ( s-n ) — 1)] and formally 



The T = random-field Ising model on a Bethe lattice with large coordination number: hysteresis and metasi 

expand W 2 (z(s, g)) around — A{g)e~ A ^ 9 ' in Eq. (1521) . which again allows us to perform 
the sum over a. Expanding the result in powers of n, we get 



2 



yd) In 2 + g[ l\z — \ z=z ( 0) 

" z< 0) p cPW~ 



EV^'^EHrOlnco^) (63) 

p=l F q=l 



where zf^ = A^e A( . The first few terms explicitly read 

4(0)2 

T^\g) = A^ 2 In 2 + 2A^A^g + r _^[ln cosh(^) - g] 

,4(0)2^(0)2 _ A (o) _ -n 

-[21ncosh(#) - lncosh(2#)] 



(1 - A(°)) 3 

1 y4(°) 3 (2/4( ) 3 - 5/4 + 6) 
+- (i _ ,4(o))5 [31ncosh(5() - 3 lncosh(2si) + lncosh(3#)]... 

(64) 

This of course gives back the power series (160]) when expanding in powers of g. To find 
an integral representation of fl63l) . we then use the identity [T5] 

, / n 2n f°° , sin(2gvrx) 

tanh(g^) = — / dx (65) 

£ J sinh(vr 2 x/5f) 

which yields by integration 

lncosh(gg) = / dx ^ X } (66) 

Jo x smh{Ti 2 x/ g) 

(the study is now restricted to g > but we know that in the end £q(<7, i?) and m*(g, H) 
are even and odd functions of g, respectively). As will be explained just below, it is 
convenient to split the domain of integration in Eq. ( 1661) into the successive intervals 
[0, 1], [1, 2], [2, 3], etc., and to change x to 1 — x in the interval [1/2, 1]. This yields 

i-l/2 

lncosh(q'g) = / /(x, g)[l — cos(2g7rx)] (67) 

Jo 

where 

/(^ = ( x + fc) sinh(7r 2 (x + £;)/#) ' (68) 

Using J2q=i(~^) q (q) = ~ 1) we then have 

V(-l)«(;)lBC08h( W ) = - / V2 rfx /(x,<?)K(l - e " 2i ™7 , (69) 

where 9ft denotes the real part, so that the series (1631) admits the following integral 
representation 

4 (0)2 

T^) = A^ 2 In 2 + 2A(°)AW 5 - 

rl/2 

+ / c/x/(x,(7)[A(°) 2 -3?{lU 2 (-AWe- A(0) e- 2 ^)}] , (70) 
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where we have interchanged summation and integration. (A priori this is not justified 



since the series Yl^Li(~ 1) P ~^^EF~ \ z -z {0) ^ — e 2l7TX ) p is n °t uniformly convergent 



for x in [0,1/2]; on the other hand, we have checked numerically that the power 
series (1611 is asymptotic to the result of Eq. (1701 for sufficiently small g, even for 
> 1.) Note that it is essential that the path in the complex plane does not 
cross the branch cuts of Wo and W-i when integrating over x in Eq. ( 1701) . The 
branch cut of Wq is {z : — oo < z < — 1/e} whereas W-\ has the double branch- 
cut {z : -oo < z < -1/e} and {z : -oo < z < 0}[Hj. Since -A^e~ A{0) > -1/e, the 
restriction of the domain of integration to the interval [0, 1/2] ensures that the imaginary 
part S{— e~ A(0) e~ 2tnx )} > 0, which explains that we use the integral representation 
(ETJ) instead of (f65]). 

Replacing T^\g) in Eq. (HHl) and f|49|) by its expression (70), and using Eqs. fl55l) 
and fl56|) . we finally obtain 

1 \i» ' / U -| o IT)* v U / 1 r> 7T)* o 



and 



where 



2 JP* v u ' 1 - 2 JP* 2 



(71) 



S Qjl (<7, if) = 2J 2 P* 2 In 2 + -(Cfo if) - g ^ 7 ) (72) 



04(0)2 

C( y , F) T«(<?) - A^ 2 In 2 - 2^^ + 

= / 1/2 rfx /(i, 9 )[i (0)2 - ^{W 2 (-A^e- AW e- 2mx )}} . (73) 
Jo 



Especially important are the limits g — > ±oo, as discussed in Ref.JTU]. It turns out 
that further analytical progress can be made in this case. Indeed, one can show that 
f(x,g) ~ g/ sin 2 (7rx) — 21n(2) when g — > +oo. As a result, it is found that 

4 72-p*2 /2-p*' 

m\(g -> +oo, if) = -ml + (1 - m* 2 ) 

1 /- 1 / 2 AW 2 - ^{W\-A^e- A[0) e- 2i ™)} 

+ o / -2/ . ? (74) 

2 Jo sm (7rx) 

and, after some additional manipulations, 

E Q>1 (g +oc,H) = ln2 dx M{W 2 (-A^e- AW e~ 2 -)} = + ^} (75) 

Jo In 6 I 

where u; = W{+A^e~ m ). 

With these expressions, we can now study the behavior of the quenched complexity 
in the field-magnetization plane and compare with the hysteresis loop. (In Appendix 
C, we indicate a change of variable in Eqs. (!73|) and (174")) which is better suited for 
numerical calculations.) 
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5. Results and discussion 



The exact equations that describes the saturation hysteresis loop on a Bethe lattice 
were obtained in Ref.[8] and the expansion of these equations to 0(1/2;) is worked out 
in Appendix B. Let us recall again that the main feature is the existence of an out- 
of-equilibrium phase transition for z > 4 in the thermodynamic limit. Whereas the 
hysteresis loop is smooth for A > A c (z), it has a jump for A < A c (z) that results 
from the flip of a finite fraction of spins (a macroscopic "avalanche"). The critical 
disorder A c (z) increases with z (see e.g. Fig. 1 in Ref.[16j) and the mean-field behavior 
described by Eq. (1421) is recovered when z — *> oo[l6]. From a mathematical point of view 
the transition for A < A c (z) is due to a saddle-node bifurcation [8]: the self-consistent 
equation (1B.2j) (that corresponds to an increasing applied field) has three real roots in 
the range ifi(A) < H < ^(A) and two of the solutions coalesce and become complex 
at the branching fields Hi and H 2 (see Fig. 3 below). As H increases, the quantity U (H) 
defined by Eq. ( IB. 21) follows initially the lowest solution and then jumps to the highest 
one at H = H 2 {A). The symmetric behavior is observed when decreasing the field. 
Whether or not the third "unstable" root and the corresponding reentrant branches 
of the hysteresis loop have a physical meaning will be discussed below. Note that an 
intermediate branch is already present in the mean-field equation (j4"2"|) below the critical 
disorder A° = Jv^fZj. Since dm^/dH < 0, this branch is also generally described 
as "unstable" , which is actually a misleading terminology since some metastable states 
may be present, as shown in Ref.pl)]. The important point for the following discussion 
is that the key quantity = 2JV(H + Jm^H)) is smaller than 1 for A > A° and on 
the stable portions of the mean-field curve for A < A° c . On the other hand, it is larger 
than 1 along the intermediate "unstable" branch (A^ = 1 at the spinodal endpoints). 



5.1. A > A° 

We first consider the strong-disorder regime in which the loop is smooth in the 
thermodynamic limit. For simplicity, we assume that A > A° > A c (z) so that the 
mean-field magnetization curve is also smooth. Then A^ < 1 and W = Wo in the 
equations of the preceding section. To illustrate numerically the general behavior, we 
take A = 1 and z = 30. 

In Fig. 1, we first show the hysteresis loop computed from the equations of Ref. 
[8] (see also Eqs. fIB.ll) and flB.2j) ) and the mean-field magnetization curve m (i?). We 
also plot the loop resulting from the expansion in 1/z, m(H) = itlq{H) + (l/z)mi(H) 
where mi(H) is given by Eqs. (1B.19j) and (1B.20|) . As can be seen, the exact result is 



accurately described by the first two terms in the 1/z expansion for the chosen value 
z = 30. 

The curves m*(g,H) = itlq{H) + (l/z)ml(g, H) versus g and Sg(m, H) = 
(I/z^q^i (m, H) versus m following the results of the preceding section are plotted 
in Figs. 2 (a) and 2 (b), respectively, for H = 0. The results are qualitatively similar for 
any other value of the field. As g varies from — oo to +oo, the magnetization m*(g, H) 
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B o- 




Figure 1. Hysteresis loop for A = 1 and z = 30 (solid black line). At the scale of 
the figure, the X/z expansion at first order (red dashed line) reproduces accurately the 
exact result. The solid blue line is the mean- field magnetization curve and the dashed 
black line is the locus of the maximum of the complexity. 




Figure 2. Magnetization m*(g,H) = m* (H) + (l/z)m*(g, H) (left panel) and 
quenched complexity Eq(to, H) = (l/z)Eg i i(m, H) (right panel) for A = 1, z = 30, 
and H = 0. The arrows in (b) indicate that g increases from — oo and +oo along the 
curve. £q(<?, H) — > when g — > ±oo so that Sq(to, H) = on the hysteresis loop (see 
text). 



increases monotonously, and the complexity, as a function of m, has a shape similar 
to that observed for z = 2|3]. In particular, the slope is infinite when g — ► ±oo, as a 
consequence of the Legendre relation g = —dY>Q(m,H)/dm. 

By studying the behavior as g — ► ±oo, we can now answer the question that 
motivated the present study: is the saturation hysteresis loop the boundary of the 
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region in the H — m plane where the quenched complexity is positive ? In the 
strong-disorder regime, the answer to this question is positive (at least at the order 
l/z). Indeed, when W = Wo, the integral in the r.h.s. of Eq. ( 174"1) is equal to 
2A(°) 2 /(1 - A®) = 8J 2 P*7(1 - 2JVq). As a result, Eq. (O identifies with Eq. 
( IB. 201) which describes the descending branch of the hysteresis loop. (Similarly, the 
ascending branch is recovered as g — > — oo.) Moreover, when w = Wo(+A^e~ A{0) ), one 
has + ^} = in Eq. (1751) so that £q,i(<? — > +oo,if), the l/z correction to the 

complexity along the hysteresis loop, is also zero. These two results prove that the curve 
E<3(m, H) = coincides with the saturation hysteresis loop at this order. (These results 
can also be obtained directly from Eq. ( 1631) by expanding formally in powers of A<-°\ 
and taking the limit g — ► +oo: one finds that T^(g) ~ 2A^°>A^g + 0(e~ 4fil ), which 
yields C(g, H) ~ 2A(°) 2 /(1-AW)^-A( ) 2 ln2 + 0(e- 49 ).) Recall that a zero complexity 
only means that the number of metastable states is not exponentially growing with N, 
but it may still vary like a power law N a with N. 

5.2. A < A c (z) 

We now turn to the low-disorder regime, illustrated by the case A = 0.3 and z = 100. 
The hysteresis loop obtained from the equations of Ref. [8] is shown in Fig. 3 (for clarity, 
the mean-field magnetization curve m (H) is only shown in Fig. 4 that zooms in the 
lower part of Fig. 3). 




Figure 3. Hysteresis loop for A = 0.3 and z — 100 (solid black line) showing jumps at 
the fields H2 and —Hi. The two intermediate branches corresponding to the "unstable" 
root of the self-consistent equation of Ref. [H] are included. The red dashed lines 
represent the "renormalized" curves H(m) = Ho(m) + (1/ z)H\(m) and the dashed 
black line is the locus of the maxima of the complexity. H\ and Hi are the branching 
fields associated with Eq. (|B.2|) (see text). 



The hysteresis loop has now a jump that occurs at H = H2 (resp. —H2) when 
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Figure 4. Zoom on the lower part of Fig. 3. The blue curve is the mean-field 
magnetization and the arrows indicate the points on this curve that correspond to 
A (0) = 1 (H = H ) and A (0) = 2 (H = H 3 ). The complexity vanishes along the two 
"stable" branches in the lower part of the loop (bold solid lines) but remains finite 
along the two green curves in the central part of the loop for < H < H3. There is an 
exponential number of metastable states in the strip between the two green curves and 
this strip is narrower than the region delimited by the dashed red lines. A first-order 
reversal curve is also shown in the lower part of the loop (solid red line) . 



increasing (resp. decreasing) the applied field. For completeness we also plot in the 
figure the reentrant branches that correspond to the "unstable" root of the self-consistent 
equation of Ref. [Ej: for Hi < H < H2 or —H% < H < —Hi when increasing or 
decreasing the field, respectively. The red dashed lines represent the "renormalized" 
curves H(m) = Ho(m) + H\(m)j z + where Ho(m) is the inverse function of itlq{H) 
(i.e. H (m (H)) = H) and Hi{m) is given by Eq. ( IB. 211) . As explained at the end 
of Appendix B, this "renormalization" of the 1/z expansion is necessary because the 
expansion m(H) = mo(H) + mi(H)/z + ... diverges when H approaches the mean-field 
branching fields, that is when — > 1. This problem can be circumvented by working 
at constant magnetization instead of constant field and expanding H instead of m. 
At the scale of Fig. 3, the resulting curves are indistinguishable from the exact ones, 
including the reentrant parts. We also display in the figure the locus of the maxima of 
the complexity which will be commented on later. 

We now discuss the results for the magnetization m*(g,H) and the quenched 
complexity Sg(m, H) = (l/^XQ^m, H). For clarity, we only focus on the lower part 
of the loop (i.e. for m < 0) as the upper one can be obtained via the symmetry 
m(—H) = —m(H). A priori, three values of the external field are noteworthy, as 
indicated in Fig. 4: the two branching fields —Hi ~ 0.391 and H2 ~ 0.413, and the 
mean-field branching field H w 0.418 (note that H 2 is close to, but smaller than, H Q ). 
As will be seen below, the field H 3 w 0.323 such that A^(H 3 ) = 2, with m* (H 3 ) being 
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Figure 5. Same as Fig. 2 for A = 0.3, z = 100, and H = 0.3 (lower part of the 
hysteresis loop). 

the "unstable" solution of Eq. (1421) . will also play a role. 

We first consider the lower part of the loop for an applied field H < —Hi. Vq is 
then computed with the solution of Eq. ( j42l) that corresponds to the lowest value of the 
magnetization mj. Then dm^jdH > and < 1 so that W = Wq in the equations of 
the preceding section. As a result, the qualitative behavior of m*(g,H) and Sg(m, H) 
is the same as in the strong-disorder regime: the region between the two branches of 
the hysteresis contains an exponential number of metastable states and the quenched 
complexity vanishes exactly along the two branches. This is illustrated in Fig. 5 (a) 
and 5 (b) for H = 0.3 (mj « —0.975 and « 0.210). It must be emphasized that 
this part of the descending branch of the hysteresis (reached in the limit g — > +oo) is 
physically accessible via a first-order reversal curve starting from the ascending branch, 
i.e. a protocol where H is first increased to a value less than H 2 and then decreased. 
Such a curve is shown in Fig. 4 (red line). States that can be reached in this way are 
called if -states [17] and it thus appears that there is an exponentially large number of 
metastable states when if -states are present (most of the metastable states however are 
not ii-states). 

The behavior is different in the central part of the loop, and we need to distinguish 
the two cases H < H% and H3 < H. The first one is illustrated in Figs. 6 (a) and 
6 (b) for H = 0. VI is then calculated with m* = so that A^ w 2.660 > 1 and 
one has to take W = W-\ in Eqs. (71)-(75). As can be seen in Fig. 6 (b), the main 
difference with the case 

A {0) 

< 1 is that the complexity does not vanish when g — ► ±oo. 
(When setting w = W^i+A^e'^) in Eq. fSJ, + is no more equal to 

zero.) Moreover, the integral in the r.h.s. of Eq. ( 1741) is not equal to 2A^ 2 /(1 — A^) so 
that Eq. (174")) does not identify with Eq. (IB. 201) . As shown in Fig. 4, the actual region 
where the number of metastable states is exponentially large, bounded by the two curves 
m*(H) = mg(-ff) + (l/z)m\{g —>■ ±oo, H) shown in green in the figure, is narrower than 
the region delimited by the reentrant branches of the hysteresis loop (red dashed lines 
in Fig. 4). This feature is in agreement with the numerical estimations of Ref. [6] for 
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Figure 6. Same as Fig. 2 for A = 0.3, z = 100, and H = (central part of the 
hysteresis loop). Note that the complexity tends to a non-zero value (« 0.0107) when 
g — ► ±oo. There are no metastable states with |m| > 0.016. 

z = 4. Therefore, these reentrant branches, which correspond to the "unstable" root of 
the self-consistent equation of Ref . [8] , bear no relation to the quenched complexity and 
do not seem to have a physical meaning. 

As H increases from and mJ(F) decreases along the reentrant part of the mean- 
field curve, decreases and the minimum value of the complexity £q(<? — > =too, H) 
also decreases. Eventually, this quantity becomes negative for « 1.9885, which 
corresponds to a field that is slightly larger than if 3 . In fact, as soon as H > H 3 , a new 
feature appears, which is illustrated in Figs. 7 (a) and 7 (b) for H = 0.36 (m^ = —0.627 
and A<® w 1.787): the magnetization m*(g,H) is no more a monotonously increasing 
function of g; in other words, the mapping {m*(g, H),T>Q(g, H)} i— > Sg(m, if) is not 
unique. This induces the awkward behavior of Sg(m, if) observed in Fig. 7 (b). 
Moreover, there is a whole range of m where Sg(m, if) is negative. 

0.015 
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X 

B 

0.005 ^"O 

W 



o 
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Figure 7. Same as Fig. 2 for A = 0.3, z = 100, and H ~ 0.36 (central part of the 
hysteresis loop). 
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Figure 8. Total complexity T^ ax {H) = {l/z)Y^ a *{H) for A = 0.3 and z = 100 (see 
note P33). As one follows the loop, the complexity increases from to ~ 0.0245 (black 
curve) and then decreases again to (red curve), as indicated by the arrows. 

The fact that dm\(g , H) / 'dg\ g=0 = (1/2) d 2 C(g, H)/dg 2 \ g=0 is negative for 1 < 
< 2 is already clear from Eq. (|60|) . the power series of T^(g) (for H = H 3 , 
m*(g,H) has an inflexion point at g = 0). Unfortunately, this signals that the solution 
of the saddle-point equations is unstable or inconsistent and that there is a flaw in 
the analytical calculations of the preceding section. Indeed, Sg(m, H) must be always 
smaller than Y% ax (H) = Y> Q (g = 0,H), which implies that <9 2 £ Q (#, H) / dg 2 \ m=cst . = 
d 2 A*W(g, H)/dg 2 = dm*(g,H)/dg must be positive at g = 0[18j. It must be also 
emphasized that the quenched complexity is a quantity that cannot become negative. 
Despite our effort, we have not succeeded in finding the origin of the problem. We 
have considered various possibilities, including a breaking of symmetry of the replica 
vector g or the existence of complex saddles [19], but unsuccessfully. Therefore, for 
H 3 < H < —Hi, we do not know the actual behavior of the complexity in the central 
region of the loop, and for —H\ < H < H2, we can only predict that the complexity 
vanishes on the lower branch of the loop (since this branch is obtained with < 1 
and W = Wq). The analytical behavior of the complexity in this intermediate region is 
obviously rather complicated. 

On the other hand, the value of the total quenched complexity Y>Q ax (H) = Sq(^ = 
0,H) is known in the whole H — m plane (cf. Eq. (1621) ). as well as the corresponding 
magnetization m max (H). These two quantities, properly "renormalized" in order to 
avoid a spurious divergence when A® -»• 1 [20J, are shown in Fig. 8 (T,Q ax (H)) and Figs. 
3 and 4 {m max (H)). T^ ax is multivalued for -H 4 < H < H 4 where H 4 = (H 2 - #0/2 
and it is maximum in the central part of the hysteresis loop. Therefore, the typical 
magnetization of the metastable states is actually a discontinuous function of H with a 
negative slope around H = 0. The surprising result that the magnetization of the most 
probable metastable states goes oppositely with the applied field is in full agreement 



The T = random-field Ising model on a Bethe lattice with large coordination number: hysteresis and metasi 

with the numerical results of Ref. [6 J (see Figs. 16 and 17 in this reference). 

Let us stress again that all the above results concern the behavior of the quenched 
complexity, associated with the typical number of metastable states. The annealed 
complexity is much more easily computed (along the lines of Ref.jl]) but much less 
informative. It indeed remains positive outside the saturation hysteresis loop for finite 
z, as shown in Ref. [I] , which signals the presence of "atypical" metastable states in this 
region. The two complexities only coincide along the mean-field magnetization curve in 
the limit z — > oo and are then equal to zero. 

6. Conclusion 

The main result of this paper concerns the organization of the metastable states of the 
RFIM in the field-magnetization plane and its relation with the saturation hysteresis 
loop. Calculations have been performed on a Bethe lattice of connectivity z in the 
large- 2 limit to order 1/z. Although some pieces of the solution are still missing, the 
following conclusions can be made, that complete and confirm earlier computations: 

1) When the hysteresis loop is smooth in the thermodynamic limit (strong-disorder 
regime), the quenched complexity Eg(m, H) is positive everywhere inside the loop, i.e. 
the number of the metastable states with a given magnetization increases exponentially 
with the system size, and it vanishes exactly along the loop. This is also the behavior 
computed analytically in one dimension [4] and observed numerically for z = 4 and on 
the cubic lattice [6]. It is therefore quite reasonable to conclude that this is a general 
result. It would be nice, of course, to have another, more direct derivation of this 
property. 

2) When the hysteresis loop is discontinuous (low-disorder regime), the quenched 
complexity is positive in two sectors: 

First, in the two regions (before and after the jump in magnetization) that can be 
reached by a field history starting from one of the saturated states. It appears that there 
is an exponentially large number of metastable states when if -states (i.e. field-reachable 
states[I7j) are present, even if most of these metastables states are not if-states. The 
complexity vanishes along the two branches of the loop that coincide with the envelope 
of the if -states. In these two regions, the typical magnetization of the metastable states 
(i.e. the magnetization of the states whose number dominates the whole distribution) 
increases with H. 

Second, in a strip of finite width (of order 1/z in our calculations) around the mean- 
field magnetization curve. This region in the middle of the loop is inaccessible to any 
field history starting from the saturated states. In addition, the complexity does not go 
continuously to zero on the borders of the strip, which moreover do not identify with 
the "unstable" branches computed from the fixed-point equations of Ref. [8]. It is in this 
strip that the number of metastable states is the largest and the typical magnetization 
decreases with H. 

The above results are in agreement with the numerical study of Ref. [6] for z = 4 
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Figure 9. Evolution of the envelope of the typical metastable states in the field- 
magnetization plane in the low-disorder regime as a function of the connectivity z: (a) 
mean-field limit (z — oo), (b) order 1/z, (c) z = 4. The complexity vanishes along the 
parts of the curves drawn in red. Certain parts are still putative (dashed lines). 

and for the cubic lattice and they suggest that in the low-disorder regime the region 
of positive (quenched) complexity evolves with connectivity as depicted in Fig. 9. We 
recall that a few metastable states are already present along the intermediate part 
of the mean-field curve [TO] . Although the precise behavior of the complexity in the 
vicinity of the "knees" is still unknown (and is certainly different on the Bethe and on 
euclidian lattices), it is now clear that the discontinuity in the hysteresis loop associated 
with a macroscopic avalanche is due to the existence of a gap in the magnetization 
of the metastable states for a range of applied field. Consequently, by controlling the 
magnetization instead of the magnetic field[9], one should observe a reentrant loop, as 
indeed observed in some magnetic materialspQ and other disordered systems [2T1 122] . 
Finally, one may wonder whether the most numerous, hence probable, states in the 
middle of the loop are accessible dynamically; it is, however, dubious that this can be 
achieved by performing a deep quench from T = ootoT = 0ata fixed field[23j. 

To conclude, we have related the out-of-equilibrium disorder-induced transition in 
the RFIM at zero temperature to the distribution of the metastable states in the field- 
magnetization plane (thus replacing the dynamic problem by a purely static calculation). 
An interesting challenge would be to describe in a similar way other out-of-equilibrium 
field-driven transitions, for instance the depinning transition [2^ [25] . In this case, one 
must take explicitely into account the geometry of the lattice. 
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Appendix A. Calculation of the complexity at the order 1/z 

In this appendix, we detail the calculation of the complexity at the order 1/z, solving 



We first note a difficulty in the 1/z expansion of the function c*(cr, r) arising from 
the presence of terms like y.r in Eq. ( J24l) . (y.r) 2 at the next order, etc... These terms 
come from the expansion of e _ *^ y T and they introduce successive derivatives of the 
Dirac 5-function when integrating over y a . This makes the integrals ill-defined (or even 
divergent) and therefore this hierarchy of equations for c*(cr,r), c^cr^r), etc... may 
just be considered as a formal and convenient way of book-keeping all the contributions 
at a certain order in 1/z. To circumvent the difficulty, however, one can simply re- 
exponentiate the problematic terms. Consider for instance the quantity J(cr, r) defined 
by Eqs. ( 1341) and ( |35l) . It can be rewritten as 



Eq. (J23D. 




a a 



1] 



(A.l) 




JV*[(l±T)l[8 K (T a ;T)-l} 



a 



2JV*Y[d K (a a ;r a ) if a ± ± 1 . 



(A.2) 



a 



where 5 k is the Kronecker symbol. This shows that the dependence of c\(ct,t) on a 
and r is rather simple. 



We now consider the quantity Xx(cr) = J2t c i( t ' <t )- F rom Eq. (|34|) . it satifies the 
self-consistent equation 




T T 

where t = Y2 a T& ' Using Eq. flA.2j) . we find after some algebra 



X 1 (±l) = ^F(r) ± JV*e-^[e 9n+Xia) - e -fl»+*i(-l)] 



T 



X x {(r) = ^F(r) + JV* e - A o[2e 9S+Xli<T) - e 9n+Xl{1) - e -fl"+*i(-l)] if a ^ ± l 



T 



(A.4) 



so that 




T 
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and 

X x {a) = Xl(1) + 2 JP*e" A S [ e 9S+x ^ - e 9n+Xl{1) ] . (A.6) 

This is the crucial self-consistency equation (145 j) . valid for all cr. 

We then calculate the quantity F(cr) defined by Eg. (1341) . After some manipulations, 
we obtain from Eqs. (21b), ([24]) and ([34} 

F{<r) = l-A{ - Xiitr) + X 2 {tr) - \x x {v)*\d&o) \{ S K (a a ; a) 

a 

+ e -K+9s+x,{a) J2 T )c*( T , CT ) 



X 



g-Ag+g.+XiCtr) /" <p( h ) dh I dxdy e iy.(x-x*) JJe(xV a )|iJAc*[l + X 1 (<r)] 

J J a 

£y a -y(i-Ac* 2 )($>°) 2 } ( A - 7 ) 



where X 2 (cr) = Er c^r, cr). The last term is computed by integration by parts over 
h, which gives 



cr 



F(cr) = [-At -X 1 (a) +X 2 (a) - ^X^af^cr) \[ 6 K (a 

a 

+ e -A5+^+^(0") ^ J( (T) T ) C *( T) CT ) 
T 

+ ae-^+s'+xWi-JAdftl + X x {a)\P* + ^(1 - Ac* 2 )V'*} JJ c^K; a) 

a 

(A.8) 

where V'(h) is the derivative of with respect to h. It thus remains to calculate 
Er / (°'> r ) c i( r ; ') and ^(±1) = Er ^(t,±1). Using Eqs. flMD and (|AT2]) . we find 
after some algebra 

^/(± l,r)c*(r,±l) = JP*[2F(±1) - ^(±1) + 2JP*e- A S ±9n+Xl{±1) ] 
r 

1(<r, r)cl(r, cr) = 2JV*F((T) + 4j2p*2 e -A 5+gs +x 1 (cr ) 

T 

= 2JV*[F(<t) - W{cr)\ ifcr#±l (A.9) 

where W(tr) = -X x (o) + X^l) - A(g) and A(g) = 2JV*e^ +gn+Xl{1 \ Inserting into 
Eq. (lA.7p . we obtain 

[1 - 2JV*e-^ ±9n+Xl{±1) ]F{± 1) = [-A* - Xi(±l) + A 2 (±l) - -X 1 (±lf]c* (±l) 
— JV*e~ A ° ±9n+Xl ^^[Xi(±l) — 2JV*e~ A * )±gn+Xl ( ± ^} 

±e -A 5±S n+X 1 (±l)|_ J Ac J[! + Xi(±l)]P* + —(1 - Ac* 2 )P*'} (A.10) 

and 

[1 - 2JP*e- A <* +ss+Xl(cr) ]F(cr) = [2Jp* e - A <* +9S+Xl(cr) ] 2 = lU 2 (cr) if a ^ ± 1 (A.ll) 
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where we have used Eq. (1A.6D and the definition of W(a). 

In order to compute X 2 (±l), we consider the normalization equation ([TO]) at the 
order 1/z 2 , which reads 

2 c o(^ t)c* 2 (t, <t) + J2 c i( ct > t K( t ' ct ) = > (A.12) 

<J T <T T 

whence 

c (l)X 2 (l)+c (-l)X 2 (-l) = -^ct(cr,T)ct(T, -) . (A.13) 

cr,r 

Using Eqs. flUD, (EOjl and flATT]) . we find 

c (l)X 2 (l) + c (-l)X 2 (-l) = F(l)[l - 2JP* e - A S+f™+^( 1 )] + F(-l)[l - 2JV*e~ A °- 9n+Xl{ - 1) ] 

-i[Xx(l) + X^-l)] 2 - -[Xx(l) + Xi(-l)][l - Jp* e -A5( e ^+X l( l) + g-sn+XU-l))] 

8 2 

_j2p*2^ e 2[-AS+gn+Xi(l)] + 5 + e -2A5+Xi(l)+Xi(-l)l 

4 T (3) (A.14) 

where we have introduced T(g) = W /2 (cr) (cf. Eq. (1521)). 

Inserting this into Eqs. flA.10l) . adding the two equations for -F(l) and F(-l), 
and using the normalization equation ( l26|) . we finally obtain Eq. ( ]44~1) . (Note that, 
remarkably, the terms involving -F(± 1) have cancelled out so that there is no need to 
compute X 2 (l) and X 2 (— 1) separately, which would imply to also consider the equation 
for c 2 (<r,r).) 

In order to calculate A^ 1 ^, one also needs the following expressions of Ac^ and 
X[ 1 \±l), solutions of Eqs. (|3"0]). (|3"T1). ([55]). and (1371) at the order n, 

A 1 ) _ — m *o 



Ar { > - a— 



2JVq] 2 



xf >(±1) = ±g(l ± mg) - . (A.15) 

Appendix B. Calculation of the hysteresis loop at the order 1/z 

In this Appendix, we compute the hysteresis loop at the order 1/z, starting from the 
equations derived by Dhar et al. [8]. The approach to the mean-field behavior has been 
studied numerically in Ref.[16j, but, as far as we know, the explicit calculation has only 
been performed in the limit z — > oo. According to Ref. [8j, the magnetization m(H) 
along the ascending branch of the hysteresis loop is given by 

l -[m{H) + 1] = QU\1 - Uy- k p k (B.l) 

fc=0 

where U(H) is solution of the self-consistent equation 

z-1 

u = Y,(r 1 )u k (i-u) z - 1 - k Pk (b.2) 

k=0 
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and Pk{H) = /J^. J / 1 _2*) V{h)dh is the probability for a spin down to flip up at the field 
H when k of its nearest neighbours are up[S] (J has been rescaled by z). In order to 
compute the 0(1/ z) correction, we assume the expansion U — Uq + t/i/z + U2/ z 2 + 
and, anticipating the fact that the sums in Eqs. fIB.ll) and flB.2p are dominated by the 
term corresponding to k — zUq, we change to the variable t = k/z and replace the sums 
by integrals that can be computed asymptotically by the Laplace method|26j. Eq. (IB. II) 
then becomes 

1 f 1 

-(m(H) + l)~z dtexp[lnQ + ztlnU + z(l -t)ln(l - U)]p(t) (B.3) 

2 jo 

with p(t) = J^°jj + jn_2t) Vifydh. Similarly, Eq. flB.21) . which is conveniently rewritten 
as U(l -U) = J2 z k=Q ^QU k (l - U) z - k p k , becomes 

U(l-U)~z f dt(l - t) exp[ln Q + ztlnU + z(l - t) ln(l - U)]p(t) . (B.4) 
Jo 

Expanding U and m(H), and using the Stirling approximation for the binomial 
coefficient in the large-z limit, we find 

' ' " ' : fl dte^f(t)[l + -g(t) + ...} 



2 (m (H) + 1) + —m\(H) + dte z ^f(t)[l + -g(t) + ...] (B.5) 

and 

U (l - U ) + -Ux(l - 2U ) + ... ~ x [± [ dte z ^\l - t)f(t)[l + -g(t) + ...] (B.6) 
where 

0(t ) =t l n ^ + ( l_ t ) ln l_^ (B.7) 

The function <p(t) is maximum for t = U , with <j)(U ) = 0, \-<p" '(U )]- 1/2 = y/U Q (l - U ) 
and f(U ) = p(U )/^yUo(l — Uq), so that, at the leading order, one obtains from Eq. 

~(to (#) + 1)=K^o) (B.10) 
and, from Eq. flB.61) . 

U Q = p(U ) (B.ll) 

Using the definition of p(t), this readily yields Eq. (|42p . the self-consistent mean-field 
equation. To compute the terms of order 1/z we need to take into account the first 
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correction to the asymptotic behavior in the Laplace'method[25]. This gives from Eq. 

das]) 



-mm = y/U {l - U ) {f(U )g(U ) - + g[<m)]2 + 2[<m)]2 

5f(U )[<f>"'(U )¥ 



24[<m)l 



(B.12) 

and, from Eq. (1B.6I) . a similar expression for U\(l— 2Uq) with /(i) replaced by (l—t)f(t). 
Using Uq = (m + l)/2, we have 



1 I — ml 



20"(^o) 8 
<f/'"(U ) l + 3m2 



?[0"(t/ o )] 2 4(1 -ml 
<1>'"{U Q ) m 



2[0"(f/ o )] 2 2 

5[<r(t/o)] 2 „ 5mg 



(B.13) 
(B.14) 
(B.15) 
(B.16) 



24[0"(f/ o )] 3 6(1 -ml) ' 
After some algebra, we obtain 

mi = AJV^U X + (1 - m 2 )J 2 V*' (B.17) 

and 

where 2JP * = 2JP(F + Jm ) = p'(c) and P *' = P'(# + Jm ), using the notations of 
sections III and IV (with m$ replaced by mo). Hence finally, 

4 72-p*2 T2<p*' 

w = -a + + (i - ^T^Jn ■ < ai9 > 

Similarly, along the descending branch (using the symmetry H — > —H, m — > — m), 

4 72-p*2 /2-p*' 

«h W = (i - "oJr^ + (i - ™o)T^jp5 • (B-20) 

For A < A° = Ja/2/7t, the function m (H) is multivalued and Eqs. (IB. 191) and 
(IB.2Q0 diverge at the mean-field spinodal where dm /dH = 2Vq/(1 — 2JVq) — > oo. 
This problem may be cured by considering the field as a function of the magnetization 
and expanding H as H(m) = H (m) + Hx(m)/z + H 2 (m)/z 2 + where H (m) 
is the inverse function of m (H), the solution of Eq. (j4"2"l) . Then Hi(m) = 
^x{H {m))/[dmo/dH\ H=Hn{m) ] = -mi(ff„(m))[l - 2JP *(#o(m)]/2P *(#o(m)) and 
from Eq. ( 1B.19I) we obtain 

1 T 2 V*'(m) 
Hx(m) = 2(1 + m) J 2 P*(m) - -(1 - m 2 ) '°l j (B.21) 



2 V ' VZ(m) 



along the ascending branch. 
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Appendix C. Alternative expression of C(g,H) 

Eqs. ( 1731) and ( 1741) are not convenient for numerical calculation because the Lambert 
function must be evaluated in the complex plane. This problem can be circumvented 
by changing the integration variable. From the equation W(z)e w ^ = z, we have 

A^e~ A(0) cos27ra; = — e u [ucosv — v sint>] 

A^e~ A ° sin27rx = e u [vcosv + usinv] (C.l) 

where u(x) and v(x) are respectively the real and imaginary parts of W{-A^e~ A{ ) e~ 2inx ). 
For W = Wq (resp. W = W^i), u is a continuous function of x that monotonically in- 
creases (resp. decreases) as x increases from to 1 and v can be univoquely expressed 
as a function of u through the relation 

v(u) = ± ( A ^ 2 e- 2{Am+ ^ - u 2 ) 1 ' 2 (C.2) 

where the signs + and — refer to W — Wq and W = W-i, respectively. A careful 
analysis shows that the inverse function x{u) is given by 

1 utsnav + v (m 

x[u) = — arctan , — A y ' < u < u\ (resp. u\ < u < —A y ') 

27r v tanf — u 

,,11 u tan v + v , , 

x(u) = — | arctan , u\ < u < Uq (resp. Uq < u < u\) 

2 2i v tant> — u 

(C.3) 

where u = u(x = 1/2) = ^{W{A ( ^e- Ai0) )} and Ul = u(x = 1/4) = ^{W{iA^ e - A{0) )}. 
Making the change of variable x — > u in Eq. (1731) . we then obtain 



C {9,H) = — du [f(x{u),g) + f{l-x(u),g] , 

27T J_ A (0) V 

(C.4) 

where the proper expression of x(u) must be used in each interval [— A^\ Ui] and 
)ui,u ) (resp. [wo, «i] and ]uu—A®]). Similarly, in Eq. (J7U), by using — dx = 
~mix C0 ^{' KX ) dx = — -jj^ cot(nx(u)) du and integrating by parts, we get 

C( 9 ,H) 2 r du 2u + * + * (C.5) 



g 7r J_a(o) tan(7rx(u)) 

when g — > +oo. Note that the integrand in Eq. (1C.4I) diverges for both u — > — A^ and 
u — > Uq when W = Wq (in both cases v — > 0), whereas it only diverges for u — » —A^> 
when W = W-\. However, these are inverse square-root singularities and the integral 
is finite, as it must be. 
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